import { distance2 } from './math';

function lerpPoint(p1, p2, t) {
	return {
		x: p1.x + (p2.x - p1.x) * t,
		y: p1.y + (p2.y - p1.y) * t
	};
}

function getPointOnCircle(x, y, cx, cy, cr) {
	const theta = Math.atan2(y - cy, x - cx);
	return {
		x: cx + Math.cos(theta) * cr,
		y: cy + Math.sin(theta) * cr
	};
}

function getPointsOnCurve(p1, cp1, cp2, p2, tolerance = 0.15) {
	// Find points on Curve (Bezier3)
	const curvePoints = [p1, cp1, cp2, p2];

	function flatness(points) {
		const p1 = points[0];
		const p2 = points[1];
		const p3 = points[2];
		const p4 = points[3];
		let ux = 3 * p2.x - 2 * p1.x - p4.x;
		ux *= ux;
		let uy = 3 * p2.y - 2 * p1.y - p4.y;
		uy *= uy;
		let vx = 3 * p3.x - 2 * p4.x - p1.x;
		vx *= vx;
		let vy = 3 * p3.y - 2 * p4.y - p1.y;
		vy *= vy;
		if (ux < vx) {
			ux = vx;
		}
		if (uy < vy) {
			uy = vy;
		}
		return ux + uy;
	}

	function splitPoints(points, tolerance, newPoints) {
		const pts = newPoints || [];
		if (flatness(points) < tolerance) {
			const p0 = points[0];
			if (pts.length) {
				const d = distance2(pts[pts.length - 1], p0);
				if (d > 1) {
					pts.push(p0);
				}
			} else {
				pts.push(p0);
			}
			pts.push(points[3]);
		} else {
			// subdivide
			const t = 0.5;
			const p1 = points[0];
			const p2 = points[1];
			const p3 = points[2];
			const p4 = points[3];
			const q1 = lerpPoint(p1, p2, t);
			const q2 = lerpPoint(p2, p3, t);
			const q3 = lerpPoint(p3, p4, t);
			const r1 = lerpPoint(q1, q2, t);
			const r2 = lerpPoint(q2, q3, t);
			const red = lerpPoint(r1, r2, t);
			splitPoints([p1, q1, r1, red], tolerance, pts);
			splitPoints([red, r2, q3, p4], tolerance, pts);
		}
		return pts;
	}

	// calc
	const results = [];
	splitPoints(curvePoints, tolerance, results);
	return results;
}

function getBezierControlPoints(keyPoints) {
	function calcCoords(arr) {
		const cp1Arr = [];
		const cp2Arr = [];
		const n = arr.length - 1;
		const a = [];
		const b = [];
		const c = [];
		const r = [];
		let m;
		let i;

		// left most segment
		a[0] = 0;
		b[0] = 2;
		c[0] = 1;
		r[0] = arr[0] + 2 * arr[1];

		// internal segments
		for (i = 1; i < n - 1; i++) {
			a[i] = 1;
			b[i] = 4;
			c[i] = 1;
			r[i] = 4 * arr[i] + 2 * arr[i + 1];
		}

		// right segment
		a[n - 1] = 2;
		b[n - 1] = 7;
		c[n - 1] = 0;
		r[n - 1] = 8 * arr[n - 1] + arr[n];

		// solves Ax=b with the Thomas algorithm (from Wikipedia)
		for (i = 1; i < n; i++) {
			m = a[i] / b[i - 1];
			b[i] = b[i] - m * c[i - 1];
			r[i] = r[i] - m * r[i - 1];
		}

		cp1Arr[n - 1] = r[n - 1] / b[n - 1];
		for (i = n - 2; i >= 0; --i) {
			cp1Arr[i] = (r[i] - c[i] * cp1Arr[i + 1]) / b[i];
		}

		// compute p2
		for (i = 0; i < n - 1; i++) {
			cp2Arr[i] = 2 * arr[i + 1] - cp1Arr[i + 1];
		}

		cp2Arr[n - 1] = 0.5 * (arr[n] + cp1Arr[n - 1]);
		return {
			cp1Arr,
			cp2Arr
		};
	}

	const xArr = [];
	const yArr = [];
	for (let i = 0; i < keyPoints.length; i++) {
		const kp = keyPoints[i];
		xArr.push(kp.x);
		yArr.push(kp.y);
	}
	const xInfo = calcCoords(xArr);
	const yInfo = calcCoords(yArr);
	const cp1Points = [];
	const cp2Points = [];
	for (let i = 0; i < xInfo.cp1Arr.length; i++) {
		cp1Points.push({
			x: xInfo.cp1Arr[i],
			y: yInfo.cp1Arr[i]
		});
		cp2Points.push({
			x: xInfo.cp2Arr[i],
			y: yInfo.cp2Arr[i]
		});
	}

	return {
		cp1Points,
		cp2Points
	};
}

function getBezierJoints(ps, t) {
	/*
	// draw bezier
	const joints = [];
	const jointStep = 0.02;
	for (let t = 0; t <= 1; t += jointStep) {
		const p = getBezierJoints(keyPoints, t);
		joints.push(p);
	}
	joints.push(keyPoints[keyPoints.length - 1]);
	*/

	let p;
	if (ps.length == 1) return ps[0];
	if (ps.length == 2) {
		p = {};
		p.x = (ps[1].x - ps[0].x) * t + ps[0].x;
		p.y = (ps[1].y - ps[0].y) * t + ps[0].y;
		return p;
	}
	if (ps.length > 2) {
		const nps = [];
		for (let i = 0; i < ps.length - 1; i++) {
			p = getBezierJoints([ps[i], ps[i + 1]], t);
			if (p) nps.push(p);
		}
		return getBezierJoints(nps, t);
	}
}

///////////////////////////////////
///////////////////////////////////
// Intersection

// https://github.com/Pomax/BezierInfo-2
// https://pomax.github.io/bezierinfo/#abc
// https://www.wolfram.com/mathematica/new-in-10/basic-and-formula-regions/curve-intersection.html
// https://github.com/thelonious/js-intersections
// https://github.com/thelonious/kld-intersections
// implemented and refiend by Du Haihang

const TOLERANCE = 1e-6;
const ACCURACY = 6;
const EPSILON = 1e-6;

function getPolyCoefs() {
	const args = arguments;
	const coefs = [];
	for (let i = args.length - 1; i >= 0; i--) coefs.push(args[i]);
	return coefs;
}

function bezout(e1, e2) {
	const AB = e1[0] * e2[1] - e2[0] * e1[1];
	const AC = e1[0] * e2[2] - e2[0] * e1[2];
	const AD = e1[0] * e2[3] - e2[0] * e1[3];
	const AE = e1[0] * e2[4] - e2[0] * e1[4];
	const AF = e1[0] * e2[5] - e2[0] * e1[5];
	const BC = e1[1] * e2[2] - e2[1] * e1[2];
	const BE = e1[1] * e2[4] - e2[1] * e1[4];
	const BF = e1[1] * e2[5] - e2[1] * e1[5];
	const CD = e1[2] * e2[3] - e2[2] * e1[3];
	const DE = e1[3] * e2[4] - e2[3] * e1[4];
	const DF = e1[3] * e2[5] - e2[3] * e1[5];
	const BFpDE = BF + DE;
	const BEmCD = BE - CD;

	return getPolyCoefs(
		AB * BC - AC * AC,
		AB * BEmCD + AD * BC - 2 * AC * AE,
		AB * BFpDE + AD * BEmCD - AE * AE - 2 * AC * AF,
		AB * DF + AD * BFpDE - 2 * AE * AF,
		AD * DF - AF * AF
	);
}

function evalPoly(coefs, x) {
	if (isNaN(x)) return null; // must be a number
	let result = 0;
	for (let i = coefs.length - 1; i >= 0; i--) result = result * x + coefs[i];
	return result;
}

function bisection(coefs, min, max) {
	let minValue = evalPoly(coefs, min);
	let maxValue = evalPoly(coefs, max);
	let result;

	if (Math.abs(minValue) <= TOLERANCE) result = min;
	else if (Math.abs(maxValue) <= TOLERANCE) result = max;
	else if (minValue * maxValue <= 0) {
		const tmp1 = Math.log(max - min);
		const tmp2 = Math.LN10 * ACCURACY;
		const iters = Math.ceil((tmp1 + tmp2) / Math.LN2);

		for (let i = 0; i < iters; i++) {
			result = 0.5 * (min + max);
			const value = evalPoly(coefs, result);

			if (Math.abs(value) <= TOLERANCE) {
				break;
			}

			if (value * minValue < 0) {
				max = result;
				maxValue = value;
			} else {
				min = result;
				minValue = value;
			}
		}
	}

	return result;
}

function getDerivative(coefs) {
	const derivativeCoefs = getPolyCoefs();
	for (let i = 1; i < coefs.length; i++) {
		derivativeCoefs.push(i * coefs[i]);
	}
	return derivativeCoefs;
}

function simplify(coefs) {
	for (let i = coefs.length - 1; i >= 0; i--) {
		if (Math.abs(coefs[i]) <= TOLERANCE) coefs.pop();
		else break;
	}
	return coefs;
}

function getLinearRoot(coefs) {
	const results = [];
	const a = coefs[1];
	if (a != 0) results.push(-coefs[0] / a);

	return results;
}

function getQuadraticRoots(coefs) {
	const results = [];

	if (coefs.length == 3) {
		const a = coefs[2];
		const b = coefs[1] / a;
		const c = coefs[0] / a;
		const d = b * b - 4 * c;

		if (d > 0) {
			const e = Math.sqrt(d);

			results.push(0.5 * (-b + e));
			results.push(0.5 * (-b - e));
		} else if (d == 0) {
			// really two roots with same value, but we only return one
			results.push(0.5 * -b);
		}
	}

	return results;
}

function getCubicRoots(coefs) {
	const results = [];
	let tmp;
	if (coefs.length == 4) {
		const c3 = coefs[3];
		const c2 = coefs[2] / c3;
		const c1 = coefs[1] / c3;
		const c0 = coefs[0] / c3;

		const a = (3 * c1 - c2 * c2) / 3;
		const b = (2 * c2 * c2 * c2 - 9 * c1 * c2 + 27 * c0) / 27;
		const offset = c2 / 3;
		let discrim = (b * b) / 4 + (a * a * a) / 27;
		const halfB = b / 2;

		if (Math.abs(discrim) <= TOLERANCE) discrim = 0;

		if (discrim > 0) {
			const e = Math.sqrt(discrim);
			let root;

			tmp = -halfB + e;
			if (tmp >= 0) root = Math.pow(tmp, 1 / 3);
			else root = -Math.pow(-tmp, 1 / 3);

			tmp = -halfB - e;
			if (tmp >= 0) root += Math.pow(tmp, 1 / 3);
			else root -= Math.pow(-tmp, 1 / 3);

			results.push(root - offset);
		} else if (discrim < 0) {
			const distance = Math.sqrt(-a / 3);
			const angle = Math.atan2(Math.sqrt(-discrim), -halfB) / 3;
			const cos = Math.cos(angle);
			const sin = Math.sin(angle);
			const sqrt3 = Math.sqrt(3);

			results.push(2 * distance * cos - offset);
			results.push(-distance * (cos + sqrt3 * sin) - offset);
			results.push(-distance * (cos - sqrt3 * sin) - offset);
		} else {
			if (halfB >= 0) tmp = -Math.pow(halfB, 1 / 3);
			else tmp = Math.pow(-halfB, 1 / 3);

			results.push(2 * tmp - offset);
			// really should return next root twice, but we return only one
			results.push(-tmp - offset);
		}
	}

	return results;
}

function getQuarticRoots(coefs) {
	const results = [];
	let d, f, t1, t2;
	if (coefs.length == 5) {
		const c4 = coefs[4];
		const c3 = coefs[3] / c4;
		const c2 = coefs[2] / c4;
		const c1 = coefs[1] / c4;
		const c0 = coefs[0] / c4;

		const resolveRoots = getCubicRoots(
			getPolyCoefs(1, -c2, c3 * c1 - 4 * c0, -c3 * c3 * c0 + 4 * c2 * c0 - c1 * c1)
		);
		const y = resolveRoots[0];
		let discrim = (c3 * c3) / 4 - c2 + y;

		if (Math.abs(discrim) <= TOLERANCE) discrim = 0;

		if (discrim > 0) {
			const e = Math.sqrt(discrim);
			t1 = (3 * c3 * c3) / 4 - e * e - 2 * c2;
			t2 = (4 * c3 * c2 - 8 * c1 - c3 * c3 * c3) / (4 * e);
			let plus = t1 + t2;
			let minus = t1 - t2;

			if (Math.abs(plus) <= TOLERANCE) plus = 0;
			if (Math.abs(minus) <= TOLERANCE) minus = 0;

			if (plus >= 0) {
				f = Math.sqrt(plus);

				results.push(-c3 / 4 + (e + f) / 2);
				results.push(-c3 / 4 + (e - f) / 2);
			}
			if (minus >= 0) {
				f = Math.sqrt(minus);

				results.push(-c3 / 4 + (f - e) / 2);
				results.push(-c3 / 4 - (f + e) / 2);
			}
		} else if (discrim < 0) {
			// no roots
		} else {
			t2 = y * y - 4 * c0;

			if (t2 >= -TOLERANCE) {
				if (t2 < 0) t2 = 0;

				t2 = 2 * Math.sqrt(t2);
				t1 = (3 * c3 * c3) / 4 - 2 * c2;
				if (t1 + t2 >= TOLERANCE) {
					d = Math.sqrt(t1 + t2);

					results.push(-c3 / 4 + d / 2);
					results.push(-c3 / 4 - d / 2);
				}
				if (t1 - t2 >= TOLERANCE) {
					d = Math.sqrt(t1 - t2);

					results.push(-c3 / 4 + d / 2);
					results.push(-c3 / 4 - d / 2);
				}
			}
		}
	}

	return results;
}

function getRoots(arr) {
	let result;
	const coefs = simplify(arr);
	switch (coefs.length - 1) {
		case 0:
			result = [];
			break;
		case 1:
			result = getLinearRoot(coefs);
			break;
		case 2:
			result = getQuadraticRoots(coefs);
			break;
		case 3:
			result = getCubicRoots(coefs);
			break;
		case 4:
			result = getQuarticRoots(coefs);
			break;
		default:
			result = [];
	}

	return result;
}

function getRootsInInterval(arr, min, max) {
	const coefs = simplify(arr);
	const roots = [];
	let root, i;

	if (coefs.length - 1 == 1) {
		root = bisection(coefs, min, max);
		if (root != null) roots.push(root);
	} else {
		// get roots of derivative
		const derivativeCoefs = getDerivative(coefs);
		const droots = getRootsInInterval(derivativeCoefs, min, max);

		if (droots.length > 0) {
			root = bisection(coefs, min, droots[0]);
			if (root != null) roots.push(root);

			for (i = 0; i <= droots.length - 2; i++) {
				root = bisection(coefs, droots[i], droots[i + 1]);
				if (root != null) roots.push(root);
			}

			root = bisection(coefs, droots[droots.length - 1], max);
			if (root != null) roots.push(root);
		} else {
			root = bisection(coefs, min, max);
			if (root != null) roots.push(root);
		}
	}
	return roots;
}

// LineLine
function computeLineLine(x1, y1, x2, y2, x3, y3, x4, y4) {
	// x1, y1, x2, y2, [line1, p1, p2]
	// x3, y3, x4, y4, [line2, p1, p2]

	// only one intersected point
	if ((x1 === x2 && y1 === y2) || (x3 === x4 && y3 === y4)) {
		return [];
	}

	const denominator = (y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1);
	if (denominator === 0) {
		// Lines are parallel
		return [];
	}

	const ua = ((x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3)) / denominator;
	const ub = ((x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3)) / denominator;

	// is the intersection along the segments
	if (ua < 0 || ua > 1 || ub < 0 || ub > 1) {
		return [];
	}

	// Return a object with the x and y coordinates of the intersection
	const x = x1 + ua * (x2 - x1);
	const y = y1 + ua * (y2 - y1);

	return [{ x, y }];
}

// LineLines - Line to Multi Polygon Lines
function computeLineLines(x1, y1, x2, y2, points) {
	let results = [];
	const len = points.length;

	for (let i = 0; i < len; i++) {
		const p1 = points[i];
		const p2 = points[(i + 1) % len];
		const inter = computeLineLine(x1, y1, x2, y2, p1.x, p1.y, p2.x, p2.y);
		results = [...results, ...inter];
	}

	return results;
}

// LinesLines - Multi Polygon Lines
function computeLinesLines(points1, points2) {
	let results = [];
	const len = points1.length;

	for (let i = 0; i < len; i++) {
		const p1 = points1[i];
		const p2 = points1[(i + 1) % len];
		const inter = computeLineLines(p1.x, p1.y, p2.x, p2.y, points2);
		results = [...results, ...inter];
	}

	return results;
}

// LineCircle
function computeLineCircle(x1, y1, x2, y2, cx, cy, cr) {
	// x1, y1, x2, y2, [line1, p1, p2]
	// cx, cy, cr, [circle center x, y, radius]
	const results = [];
	const v1x = x2 - x1;
	const v1y = y2 - y1;
	const v2x = x1 - cx;
	const v2y = y1 - cy;

	let b;
	let c;
	let d;
	let u1;
	let u2;
	b = v1x * v2x + v1y * v2y;
	c = 2 * (v1x * v1x + v1y * v1y);
	b *= -2;
	d = Math.sqrt(b * b - 2 * c * (v2x * v2x + v2y * v2y - cr * cr));
	if (isNaN(d)) {
		return results; // no intersection
	}

	// the unit distance of point one and two on the line
	u1 = (b - d) / c;
	u2 = (b + d) / c;
	if (u1 <= 1 && u1 >= 0) {
		// add point if on the line
		results.push({
			x: x1 + v1x * u1,
			y: y1 + v1y * u1
		});
	}
	if (u2 <= 1 && u2 >= 0) {
		// second point if on the line
		results.push({
			x: x1 + v1x * u2,
			y: y1 + v1y * u2
		});
	}
	return results;
}

// LineEllipse
function computeLineEllipse(x1, y1, x2, y2, cx, cy, rx, ry) {
	const results = [];
	const l1 = {
		x: x1,
		y: y1
	};
	const l2 = {
		x: x2,
		y: y2
	};
	const dir = {
		x: x2 - x1,
		y: y2 - y1
	};
	const diff = {
		x: x1 - cx,
		y: y1 - cy
	};
	const mDir = { x: dir.x / (rx * rx), y: dir.y / (ry * ry) };
	const mDiff = { x: diff.x / (rx * rx), y: diff.y / (ry * ry) };

	const a = dir.x * mDir.x + dir.y * mDir.y;
	const b = dir.x * mDiff.x + dir.y * mDiff.y;
	const c = diff.x * mDiff.x + diff.y * mDiff.y - 1.0;
	const d = b * b - a * c;

	if (d < 0) {
		// Outside
	} else if (d > 0) {
		const root = Math.sqrt(d);
		const t_a = (-b - root) / a;
		const t_b = (-b + root) / a;

		if ((t_a < 0 || 1 < t_a) && (t_b < 0 || 1 < t_b)) {
			if ((t_a < 0 && t_b < 0) || (t_a > 1 && t_b > 1)) {
				// Outside
			} else {
				// Inside
			}
		} else {
			if (0 <= t_a && t_a <= 1) {
				results.push(lerpPoint(l1, l2, t_a));
			}
			if (0 <= t_b && t_b <= 1) {
				results.push(lerpPoint(l1, l2, t_b));
			}
		}
	} else {
		const t = -b / a;
		if (0 <= t && t <= 1) {
			results.push(lerpPoint(l1, l2, t));
		} else {
			// Outside
		}
	}

	return results;
}

// LineRect
function computeLineRect(x1, y1, x2, y2, tx, ty, tw, th) {
	const inter1 = computeLineLine(x1, y1, x2, y2, tx, ty, tx + tw, ty);
	const inter2 = computeLineLine(x1, y1, x2, y2, tx + tw, ty, tx + tw, ty + th);
	const inter3 = computeLineLine(x1, y1, x2, y2, tx, ty + th, tx + tw, ty + th);
	const inter4 = computeLineLine(x1, y1, x2, y2, tx, ty, tx, ty + th);

	let results = [];
	results = [...results, ...inter1];
	results = [...results, ...inter2];
	results = [...results, ...inter3];
	results = [...results, ...inter4];
	return results;
}

// EllipseEllipse
function computeEllipseEllipse(x1, y1, r1x, r1y, x2, y2, r2x, r2y) {
	// x1, y1, r1x, r1y, [circle1 center x, y, radiusX, radiusY]
	// x2, y2, r2x, r2y, [circle2 center x, y, radiusX, radiusY]

	const a = [
		r1y * r1y,
		0,
		r1x * r1x,
		-2 * r1y * r1y * x1,
		-2 * r1x * r1x * y1,
		r1y * r1y * x1 * x1 + r1x * r1x * y1 * y1 - r1x * r1x * r1y * r1y
	];
	const b = [
		r2y * r2y,
		0,
		r2x * r2x,
		-2 * r2y * r2y * x2,
		-2 * r2x * r2x * y2,
		r2y * r2y * x2 * x2 + r2x * r2x * y2 * y2 - r2x * r2x * r2y * r2y
	];

	const yPolyCoefs = bezout(a, b);
	const yRoots = getRoots(yPolyCoefs);
	const epsilon = 1e-3;
	const norm0 = (a[0] * a[0] + 2 * a[1] * a[1] + a[2] * a[2]) * epsilon;
	const norm1 = (b[0] * b[0] + 2 * b[1] * b[1] + b[2] * b[2]) * epsilon;
	const results = [];

	for (let y = 0; y < yRoots.length; y++) {
		const xPolyCoefs = getPolyCoefs(a[0], a[3] + yRoots[y] * a[1], a[5] + yRoots[y] * (a[4] + yRoots[y] * a[2]));
		const xRoots = getRoots(xPolyCoefs);

		for (let x = 0; x < xRoots.length; x++) {
			let test =
				(a[0] * xRoots[x] + a[1] * yRoots[y] + a[3]) * xRoots[x] + (a[2] * yRoots[y] + a[4]) * yRoots[y] + a[5];
			if (Math.abs(test) < norm0) {
				test =
					(b[0] * xRoots[x] + b[1] * yRoots[y] + b[3]) * xRoots[x] +
					(b[2] * yRoots[y] + b[4]) * yRoots[y] +
					b[5];
				if (Math.abs(test) < norm1) {
					results.push({
						x: xRoots[x],
						y: yRoots[y]
					});
				}
			}
		}
	}

	return results;
}

// CircleCircle
function computeCircleCircle(x1, y1, r1, x2, y2, r2) {
	const results = [];
	const c1 = { x: x1, y: y1 };
	const c2 = { x: x2, y: y2 };
	const r_max = r1 + r2;
	const r_min = Math.abs(r1 - r2);
	const dx = x1 - x2;
	const dy = y1 - y2;
	const c_dist = Math.sqrt(dx * dx + dy * dy);

	if (c_dist > r_max) {
		// Outside
	} else if (c_dist < r_min) {
		// Inside
	} else {
		const a = (r1 * r1 - r2 * r2 + c_dist * c_dist) / (2 * c_dist);
		const h = Math.sqrt(r1 * r1 - a * a);
		const p = lerpPoint(c1, c2, a / c_dist);
		const b = h / c_dist;

		results.push({
			x: p.x - b * (c2.y - c1.y),
			y: p.y + b * (c2.x - c1.x)
		});
		results.push({
			x: p.x + b * (c2.y - c1.y),
			y: p.y - b * (c2.x - c1.x)
		});
	}
	return results;
}

// CircleEllipse
function computeCircleEllipse(x1, y1, r1, x2, y2, r2x, r2y) {
	return computeEllipseEllipse(x1, y1, r1, r1, x2, y2, r2x, r2y);
}

// LineCurve (Bezier3)
function computeLineCurve(ax, ay, bx, by, p1x, p1y, p2x, p2y, p3x, p3y, p4x, p4y) {
	// Line, ax, ay, bx, by, [line.p1, line.p2]
	// Curve, [head, cp1, cp2, foot];

	let a; // temporary variables

	let b;
	let c;
	let c3; // coefficients of cubic
	let c2;
	let c1;
	let c0;
	let cl; // c coefficient for normal form of line
	let n; // normal for normal form of line
	const min = {
		x: Math.min(ax, bx),
		y: Math.min(ay, by)
	};
	const max = {
		x: Math.max(ax, bx),
		y: Math.max(ay, by)
	};
	const p1 = { x: p1x, y: p1y };
	const p2 = { x: p2x, y: p2y };
	const p3 = { x: p3x, y: p3y };
	const p4 = { x: p4x, y: p4y };
	const results = [];

	// Calculate the coefficients
	a = { x: p1x * -1, y: p1y * -1 };
	b = { x: p2x * 3, y: p2y * 3 };
	c = { x: p3x * -3, y: p3y * -3 };
	c3 = {
		x: a.x + b.x + c.x + p4.x,
		y: a.y + b.y + c.y + p4.y
	};

	a = { x: p1x * 3, y: p1y * 3 };
	b = { x: p2x * -6, y: p2y * -6 };
	c = { x: p3x * 3, y: p3y * 3 };
	c2 = {
		x: a.x + b.x + c.x,
		y: a.y + b.y + c.y
	};

	a = { x: p1x * -3, y: p1y * -3 };
	b = { x: p2x * 3, y: p2y * 3 };
	c1 = {
		x: a.x + b.x,
		y: a.y + b.y
	};

	c0 = p1;

	// Convert line to normal form: ax + by + c = 0
	// Find normal to line: negative inverse of original line's slope
	n = { x: ay - by, y: bx - ax };
	cl = ax * by - bx * ay;

	const roots = getRoots(
		getPolyCoefs(
			n.x * c3.x + n.y * c3.y,
			n.x * c2.x + n.y * c2.y,
			n.x * c1.x + n.y * c1.y,
			n.x * c0.x + n.y * c0.y + cl
		)
	);

	for (let i = 0; i < roots.length; i++) {
		const t = roots[i];

		if (0 <= t && t <= 1) {
			// We're within the Bezier curve
			// Find point on Bezier
			const p5 = lerpPoint(p1, p2, t);
			const p6 = lerpPoint(p2, p3, t);
			const p7 = lerpPoint(p3, p4, t);
			const p8 = lerpPoint(p5, p6, t);
			const p9 = lerpPoint(p6, p7, t);
			const p10 = lerpPoint(p8, p9, t);

			// See if point is on line segment
			// Had to make special cases for vertical and horizontal lines due
			// to slight errors in calculation of p10

			if (ax == bx) {
				if (min.y <= p10.y && p10.y <= max.y) {
					results.push(p10);
				}
			} else if (ay == by) {
				if (min.x <= p10.x && p10.x <= max.x) {
					results.push(p10);
				}
			} else if (p10.x >= min.x && p10.y >= min.y && p10.x <= max.x && p10.y <= max.y) {
				results.push(p10);
			}
		}
	}

	return results;
}

// EllipseCurve (Bezier3)
function computeEllipseCurve(cx, cy, rx, ry, p1x, p1y, p2x, p2y, p3x, p3y, p4x, p4y) {
	// Ellipse [centerx, centerx, radiusx, radiusy];
	// Curve, [head, cp1, cp2, foot];
	const results = [];

	const c3 = { x: p1x * -1 + p2x * 3 + p3x * -3 + p4x, y: p1y * -1 + p2y * 3 + p3y * -3 + p4y };
	const c2 = { x: p1x * 3 + p2x * -6 + p3x * 3, y: p1y * 3 + p2y * -6 + p3y * 3 };
	const c1 = { x: p1x * -3 + p2x * 3, y: p1y * -3 + p2y * 3 };
	const c0 = { x: p1x, y: p1y };

	const rxrx = rx * rx;
	const ryry = ry * ry;
	const polyCoefs = getPolyCoefs(
		c3.x * c3.x * ryry + c3.y * c3.y * rxrx,
		2 * (c3.x * c2.x * ryry + c3.y * c2.y * rxrx),
		2 * (c3.x * c1.x * ryry + c3.y * c1.y * rxrx) + c2.x * c2.x * ryry + c2.y * c2.y * rxrx,
		2 * c3.x * ryry * (c0.x - cx) + 2 * c3.y * rxrx * (c0.y - cy) + 2 * (c2.x * c1.x * ryry + c2.y * c1.y * rxrx),
		2 * c2.x * ryry * (c0.x - cx) + 2 * c2.y * rxrx * (c0.y - cy) + c1.x * c1.x * ryry + c1.y * c1.y * rxrx,
		2 * c1.x * ryry * (c0.x - cx) + 2 * c1.y * rxrx * (c0.y - cy),
		c0.x * c0.x * ryry -
			2 * c0.y * cy * rxrx -
			2 * c0.x * cx * ryry +
			c0.y * c0.y * rxrx +
			cx * cx * ryry +
			cy * cy * rxrx -
			rxrx * ryry
	);
	const roots = getRootsInInterval(polyCoefs, 0, 1);

	for (const t of roots) {
		const pt = { x: 0, y: 0 };
		pt.x = c3.x * t * t * t + c2.x * t * t + c1.x * t + c0.x;
		pt.y = c3.y * t * t * t + c2.y * t * t + c1.y * t + c0.y;
		results.push(pt);
	}

	return results;
}

// CurveCurve (Bezier3) faster than js-intersections polycoefs
function computeCurveCurve(a1x, a1y, a2x, a2y, a3x, a3y, a4x, a4y, p1x, p1y, p2x, p2y, p3x, p3y, p4x, p4y) {
	const results = [];
	const pts1 = [a1x, a1y, a2x, a2y, a3x, a3y, a4x, a4y];
	const pts2 = [p1x, p1y, p2x, p2y, p3x, p3y, p4x, p4y];

	let tMin = 0.0;
	let tMax = 1.0;
	let uMin = 0.0;
	let uMax = 1.0;
	let oldTDiff = 1.0;
	let reverse = false;
	let recursion = 0;
	let recursionLimit = 32;
	let tLimit = 0.8;

	function isZero(val) {
		return Math.abs(val) <= EPSILON;
	}

	function signedDistance(px, py, vx, vy, x, y) {
		vx -= px;
		vy -= py;
		if (isZero(vx)) {
			return vy >= 0 ? px - x : x - px;
		} else if (isZero(vy)) {
			return vx >= 0 ? y - py : py - y;
		} else {
			return (vx * (y - py) - vy * (x - px)) / Math.sqrt(vx * vx + vy * vy);
		}
	}

	function convexHull(dq0, dq1, dq2, dq3) {
		let p0 = [0, dq0];
		let p1 = [1.0 / 3, dq1];
		let p2 = [2.0 / 3, dq2];
		let p3 = [1, dq3];

		// Find signed distance of p1 and p2 from line [ p0, p3 ]
		let dist1 = signedDistance(0, dq0, 1, dq3, 1.0 / 3, dq1);
		let dist2 = signedDistance(0, dq0, 1, dq3, 2.0 / 3, dq2);

		let flip = false;
		let hull;

		if (dist1 * dist2 < 0) {
			hull = [
				[p0, p1, p3],
				[p0, p2, p3]
			];
			flip = dist1 < 0;
		} else {
			let pmax;
			let cross = 0;
			let distZero = dist1 == 0 || dist2 == 0;
			if (Math.abs(dist1) > Math.abs(dist2)) {
				pmax = p1;
				cross = ((dq3 - dq2 - (dq3 - dq0) / 3.0) * (2 * (dq3 - dq2) - dq3 + dq1)) / 3.0;
			} else {
				pmax = p2;
				// apex is dq0 in this case, and the other apex point is dq3
				// vector dqapex -> dqapex2 or base vector which is already part
				// of the hull.
				cross = ((dq1 - dq0 + (dq0 - dq3) / 3.0) * (-2 * (dq0 - dq1) + dq0 - dq2)) / 3.0;
			}

			// Compare cross products of these vectors to determine if the point
			// is in the triangle [ p3, pmax, p0 ], or if it is a quadrilateral.
			hull =
				cross < 0 || distZero
					? [
							[p0, pmax, p3],
							[p0, p3]
					  ]
					: [
							[p0, p1, p2, p3],
							[p0, p3]
					  ];
			flip = dist1 ? dist1 < 0 : dist2 < 0;
		}
		if (flip) {
			hull.reverse();
		}
		return hull;
	}

	function clipConvexHull(hullTop, hullBottom, dMin, dMax) {
		if (hullTop[0][1] < dMin) {
			return clipConvexHullPart(hullTop, true, dMin);
		} else if (hullBottom[0][1] > dMax) {
			return clipConvexHullPart(hullBottom, false, dMax);
		} else {
			return hullTop[0][0];
		}
	}

	function clipConvexHullPart(part, top, threshold) {
		let [px, py] = part[0];
		for (let i = 1; i < part.length; i++) {
			let [qx, qy] = part[i];
			if (top ? qy >= threshold : qy <= threshold) {
				return px + ((threshold - py) * (qx - px)) / (qy - py);
			}
			[px, py] = [qx, qy];
		}
		// All points of hull are above / below the threshold
		return null;
	}

	function getFatline(v) {
		// Starting point of the curve
		let q0x = v[0];
		let q0y = v[1];
		// End point of the curve
		let q3x = v[6];
		let q3y = v[7];
		// Calculate the fat-line L, for Q is the baseline l and two
		// offsets which completely encloses the curve P.
		let d1 = signedDistance(q0x, q0y, q3x, q3y, v[2], v[3]) || 0;
		let d2 = signedDistance(q0x, q0y, q3x, q3y, v[4], v[5]) || 0;
		let factor = d1 * d2 > 0 ? 3.0 / 4.0 : 4.0 / 9.0; // Get a tighter fit
		let dMin = factor * Math.min(0, d1, d2);
		let dMax = factor * Math.max(0, d1, d2);
		// The width of the 'fatline' is |dMin| + |dMax|
		return [dMin, dMax];
	}

	function subdivide(v, t = 0.5) {
		let [p1x, p1y, c1x, c1y, c2x, c2y, p2x, p2y] = v;

		// Triangle computation, with loops unrolled.
		let u = 1 - t;
		// Interpolate from 4 to 3 points
		let p3x = u * p1x + t * c1x;
		let p3y = u * p1y + t * c1y;
		let p4x = u * c1x + t * c2x;
		let p4y = u * c1y + t * c2y;
		let p5x = u * c2x + t * p2x;
		let p5y = u * c2y + t * p2y;
		// Interpolate from 3 to 2 points
		let p6x = u * p3x + t * p4x;
		let p6y = u * p3y + t * p4y;
		let p7x = u * p4x + t * p5x;
		let p7y = u * p4y + t * p5y;
		// Interpolate from 2 points to 1 point
		let p8x = u * p6x + t * p7x;
		let p8y = u * p6y + t * p7y;

		// We now have all the values we need to build the sub-curves [left, right]:
		return [
			[p1x, p1y, p3x, p3y, p6x, p6y, p8x, p8y],
			[p8x, p8y, p7x, p7y, p5x, p5y, p2x, p2y]
		];
	}

	function getPart(v, t1, t2) {
		if (t1 > 0) {
			v = subdivide(v, t1)[1]; // right
		}
		// Interpolate the parameter at 't2' in the new curve and cut there.
		if (t2 < 1) {
			v = subdivide(v, (t2 - t1) / (1.0 - t1))[0]; // left
		}
		return v;
	}

	function evaluate(v, t, type) {
		let [p1x, p1y, c1x, c1y, c2x, c2y, p2x, p2y] = v;
		let x, y;

		// Handle special case at beginning / end of curve
		if (type === 0 && (t < TOLERANCE || t > 1 - TOLERANCE)) {
			let isZero = t < TOLERANCE;
			x = isZero ? p1x : p2x;
			y = isZero ? p1y : p2y;
		} else {
			// Calculate the polynomial coefficients.
			let cx = 3.0 * (c1x - p1x);
			let bx = 3.0 * (c2x - c1x) - cx;
			let ax = p2x - p1x - cx - bx;

			let cy = 3.0 * (c1y - p1y);
			let by = 3.0 * (c2y - c1y) - cy;
			let ay = p2y - p1y - cy - by;

			if (type === 0) {
				// Calculate the curve point at parameter value t
				x = ((ax * t + bx) * t + cx) * t + p1x;
				y = ((ay * t + by) * t + cy) * t + p1y;
			} else {
				// 1: tangent, 1st derivative
				// 2: normal, 1st derivative
				// 3: curvature, 1st derivative & 2nd derivative
				// Prevent tangents and normals of length 0:
				// http:#stackoverflow.com/questions/10506868/
				if ((t < TOLERANCE && c1x == p1x && c1y == p1y) || (t > 1 - TOLERANCE && c2x == p2x && c2y == p2y)) {
					x = c2x - c1x;
					y = c2y - c1y;
				} else if (t < TOLERANCE) {
					x = cx;
					y = cy;
				} else if (t > 1 - TOLERANCE) {
					x = 3.0 * (p2x - c2x);
					y = 3.0 * (p2y - c2y);
				} else {
					// Simply use the derivation of the bezier function for both
					// the x and y coordinates:
					x = (3.0 * ax * t + 2 * bx) * t + cx;
					y = (3.0 * ay * t + 2 * by) * t + cy;
				}
				if (type === 3) {
					// Calculate 2nd derivative, and curvature from there:
					// http://cagd.cs.byu.edu/~557/text/ch2.pdf page#31
					// k = |dx * d2y - dy * d2x| / (( dx^2 + dy^2 )^(3/2))
					let x2 = 6.0 * ax * t + 2.0 * bx;
					let y2 = 6.0 * ay * t + 2.0 * by;
					return (x * y2 - y * x2) / Math.pow(x * x + y * y, 3.0 / 2);
				}
			}
		}
		// The normal is simply the rotated tangent:
		return type == 2 ? [y, -x] : [x, y];
	}

	function projectCurves(v1, v2, tMin, tMax, uMin, uMax, oldTDiff, reverse, recursion, recursionLimit, tLimit) {
		if (recursion > recursionLimit) {
			return [];
		}

		// Let P be the first curve and Q be the second
		let q0x = v2[0];
		let q0y = v2[1];
		let q3x = v2[6];
		let q3y = v2[7];

		// Calculate the fat-line L for Q is the baseline l and two
		// offsets which completely encloses the curve P.
		let [dMin, dMax] = getFatline(v2);

		// Calculate non-parametric bezier curve D(ti, di(t)) - di(t) is the
		// distance of P from the baseline l of the fat-line, ti is equally
		// spaced in [0, 1]
		let dp0 = signedDistance(q0x, q0y, q3x, q3y, v1[0], v1[1]);
		let dp1 = signedDistance(q0x, q0y, q3x, q3y, v1[2], v1[3]);
		let dp2 = signedDistance(q0x, q0y, q3x, q3y, v1[4], v1[5]);
		let dp3 = signedDistance(q0x, q0y, q3x, q3y, v1[6], v1[7]);
		let tMinNew = 0.0;
		let tMaxNew = 0.0;
		let tDiff = 0.0;

		// NOTE: the recursion threshold of 4 is needed to prevent issue #571
		// from occurring: https://github.com/paperjs/paper.js/issues/571
		if (q0x == q3x && uMax - uMin <= EPSILON && recursion > 4) {
			// The fatline of Q has converged to a point, the clipping is not
			// reliable. Return the value we have even though we will miss the
			// precision.
			tMaxNew = tMinNew = (tMax + tMin) / 2.0;
			tDiff = 0;
		} else {
			// Get the top and bottom parts of the convex-hull
			let hull = convexHull(dp0, dp1, dp2, dp3);
			let top = hull[0];
			let bottom = hull[1];
			let tMinClip;
			let tMaxClip;
			// Clip the convex-hull with dMin and dMax
			tMinClip = clipConvexHull(top, bottom, dMin, dMax);
			top.reverse();
			bottom.reverse();
			tMaxClip = clipConvexHull(top, bottom, dMin, dMax);
			// No intersections if one of the tvalues are null or 'undefined'
			if (tMinClip === null || tMaxClip === null) {
				return [];
			}
			// Clip P with the fatline for Q
			let v1d = v1;
			v1 = getPart(v1, tMinClip, tMaxClip);
			tDiff = tMaxClip - tMinClip;
			// tMin and tMax are within the range (0, 1). We need to project it
			// to the original parameter range for v2.
			tMinNew = tMax * tMinClip + tMin * (1 - tMinClip);
			tMaxNew = tMax * tMaxClip + tMin * (1 - tMaxClip);
		}

		let res;

		// Check if we need to subdivide the curves
		if (oldTDiff > tLimit && tDiff > tLimit) {
			// Subdivide the curve which has converged the least.
			if (tMaxNew - tMinNew > uMax - uMin) {
				let parts = subdivide(v1, 0.5);
				let t = tMinNew + (tMaxNew - tMinNew) / 2.0;
				res = [].concat(
					projectCurves(
						v2,
						parts[0],
						uMin,
						uMax,
						tMinNew,
						t,
						tDiff,
						!reverse,
						recursion + 1,
						recursionLimit,
						tLimit
					),
					projectCurves(
						v2,
						parts[1],
						uMin,
						uMax,
						t,
						tMaxNew,
						tDiff,
						!reverse,
						recursion + 1,
						recursionLimit,
						tLimit
					)
				);
			} else {
				let parts = subdivide(v2, 0.5);
				let t = uMin + (uMax - uMin) / 2.0;
				res = [].concat(
					projectCurves(
						parts[0],
						v1,
						uMin,
						t,
						tMinNew,
						tMaxNew,
						tDiff,
						!reverse,
						recursion + 1,
						recursionLimit,
						tLimit
					),
					projectCurves(
						parts[1],
						v1,
						t,
						uMax,
						tMinNew,
						tMaxNew,
						tDiff,
						!reverse,
						recursion + 1,
						recursionLimit,
						tLimit
					)
				);
			}
		} else if (Math.max(uMax - uMin, tMaxNew - tMinNew) < TOLERANCE) {
			// We have isolated the intersection with sufficient precision
			let t1 = tMinNew + (tMaxNew - tMinNew) / 2.0;
			let t2 = uMin + (uMax - uMin) / 2.0;
			if (reverse) {
				res = [[t2, evaluate(v2, t2, 0), t1, evaluate(v1, t1, 0)]];
			} else {
				res = [[t1, evaluate(v1, t1, 0), t2, evaluate(v2, t2, 0)]];
			}
		} else {
			res = projectCurves(
				v2,
				v1,
				uMin,
				uMax,
				tMinNew,
				tMaxNew,
				tDiff,
				!reverse,
				recursion + 1,
				recursionLimit,
				tLimit
			);
		}

		return res;
	}

	const inters = projectCurves(
		pts1,
		pts2,
		tMin,
		tMax,
		uMin,
		uMax,
		oldTDiff,
		reverse,
		recursion,
		recursionLimit,
		tLimit
	);

	for (const inter of inters) {
		if (inter) {
			const [t1, p1, t2, p2] = inter;
			if (p1) results.push({ x: p1[0], y: p1[1] });
			if (p2) results.push({ x: p2[0], y: p2[1] });
		}
	}

	return results;
}

export {
	// Methods
	lerpPoint,
	getPointOnCircle,
	getPointsOnCurve,
	getBezierControlPoints,
	getBezierJoints,
	// Intersection
	computeLineLine,
	computeLineLines,
	computeLinesLines,
	computeLineCircle,
	computeLineEllipse,
	computeLineRect,
	computeEllipseEllipse,
	computeCircleCircle,
	computeCircleEllipse,
	computeLineCurve,
	computeEllipseCurve,
	computeCurveCurve
};
